A Markov model is a powerful tool for analyzing a patient’s journey over an extended period. Unlike a decision tree, which is a snapshot in time, a Markov model simulates how a group of patients moves through different health states over many years. It is perfect for analyzing chronic diseases like heart disease, where a patient’s condition can change over time. In this exercise, we will build a Markov model to compare two treatment strategies for coronary artery disease: Percutaneous Coronary Intervention (PCI) and Coronary Artery Bypass Grafting (CABG). We will simulate a cohort of patients over a 10-year period, tracking their health states, costs, and quality-adjusted life years (QALYs). —
1.1 Key Concepts
Health States: These are the possible clinical conditions a patient can be in. We will use a simplified set of states for our model:
Stable: The patient is healthy and has stable symptoms.
Myocardial Infarction (MI): The patient has a heart attack.
Reintervention: The patient requires a repeat procedure.
Death: The patient has died (this is an absorbing state, meaning once a patient enters this state, they cannot leave).
Cycles: The model runs in discrete time periods, or cycles. Our cycles will be one year long.
Transition Probabilities: These are the probabilities that a patient will move from one health state to another in a single cycle. These values come from a literature review of clinical trials.
Costs and Utilities: Each health state has associated costs (e.g., medical expenses) and utilities (e.g., quality of life). We will use these to calculate the total cost and QALYs for each treatment strategy.
1.2 Visualizing the Markov Model
This diagram shows how a patient can move between the different health states in our model. The arrows represent the possible transitions in a single cycle.
Code
# Check and install DiagrammeR if not already installedif (!requireNamespace("DiagrammeR", quietly =TRUE)) {install.packages("DiagrammeR")}library(DiagrammeR)# Create the visual representation of the Markov modelgrViz(" digraph G { graph [overlap = true, splines = true, fontname = Helvetica, layout = dot] // Define nodes with shapes and labels node [shape = rectangle, style = filled, color = LightBlue, fontname = Helvetica] Stable; MI [color = LightCoral]; Reintervention [color = LightYellow]; Death [shape = hexagon, color = Gray]; // Define transitions (edges) Stable -> Stable [label = 'P(stable to stable)']; Stable -> MI [label = 'P(stable to MI)']; Stable -> Reintervention [label = 'P(stable to reintervention)']; Stable -> Death [label = 'P(stable to death)']; MI -> Stable [label = 'P(MI to stable)']; MI -> Reintervention [label = 'P(MI to reintervention)']; MI -> Death [label = 'P(MI to death)']; Reintervention -> Stable [label = 'P(reintervention to stable)']; Reintervention -> Death [label = 'P(reintervention to death)']; // Death is an absorbing state Death -> Death [label = 'P(death to death) = 1', color = Gray]; }")
2 Building a Markov Model Step-by-Step
Objective: To build a working Markov model in R to compare the long-term cost and health outcomes of PCI versus CABG.
2.1 Step 1: Define Model Parameters
In this first step, we’ll define all the parameters for our model. These are the costs, utilities, and transition probabilities we’ve gathered from our hypothetical literature review.
# Costs in Indian Rupees (INR)cost_pci_initial <-200000# Cost of the initial PCI procedure for the cohortcost_cabg_initial <-300000# Cost of the initial CABG procedure for the cohortcost_stable <-50000# Annual cost for a stable patientcost_mi <-20000# Cost of a myocardial infarction eventcost_reintervention <-150000# Cost of a reinterventioncost_death <-0# Cost for the absorbing state
2.1.3 1.2.1 Show costs
Code
costs_table <- tibble::tibble(State = state_names,Cost =c(cost_stable, cost_mi, cost_reintervention, cost_death))# Check and install gt if not already installed for table formattingif (!requireNamespace("gt", quietly =TRUE)) {install.packages("gt")}library(gt)gt(costs_table)
State
Cost
Stable
50000
MI
20000
Reintervention
150000
Death
0
2.1.4 1.3 Utilities
Code
# Utilities (Quality-Adjusted Life Years - QALYs)utility_stable <-0.85# Utility value for a stable patientutility_mi <-0.6# Utility value during the year of an MIutility_reintervention <-0.7# Utility value during the year of a reinterventionutility_death <-0# Utility value for the absorbing state
# Discountingdiscount_rate_cost <-0.03# A common discount rate for costs (5%)discount_rate_qaly <-0.03# A common discount rate for QALYs (5%)
2.1.7 1.5 Define Transition Probabilities for PCI and CABG
2.1.7.1 1.5.1 PCI Transition Probabilities
Code
# Transition probabilities from literature for a 1-year cyclepci_trans <-list(stable_to_stable =0.70,stable_to_mi =0.15,stable_to_reintervention =0.10,stable_to_death =0.05,mi_to_reintervention =0.20,mi_to_stable =0.50,mi_to_death =0.30,reintervention_to_stable =0.80,reintervention_to_death =0.20)
The transition matrix is the engine of our Markov model. It’s a square matrix where each row represents the starting health state and each column represents the destination health state. The value at each intersection is the probability of that transition occurring.
Instructions: Create two matrices, one for PCI and one for CABG, using the probabilities we defined in the previous step.
2.2.1 2.1 Matrix Construction for PCI
Code
# PCI Transition Matrixtrans_pci <-matrix(c( pci_trans$stable_to_stable, pci_trans$stable_to_mi, pci_trans$stable_to_reintervention, pci_trans$stable_to_death, # From Stable... pci_trans$mi_to_stable, 0, pci_trans$mi_to_reintervention, pci_trans$mi_to_death, # From MI... pci_trans$reintervention_to_stable, 0, 0, pci_trans$reintervention_to_death, # From Reintervention...0, 0, 0, 1# From Death (absorbing state) ),nrow =4,byrow =TRUE,dimnames =list(state_names, state_names))# Print the matrices to see their structureprint("PCI Transition Matrix:")
[1] "PCI Transition Matrix:"
Code
print(trans_pci)
Stable MI Reintervention Death
Stable 0.7 0.15 0.1 0.05
MI 0.5 0.00 0.2 0.30
Reintervention 0.8 0.00 0.0 0.20
Death 0.0 0.00 0.0 1.00
2.2.2 2.1.1 Validation of Transition Matrices
Code
rowSums(trans_pci)
Stable MI Reintervention Death
1 1 1 1
Code
trans_pci["Death", ]
Stable MI Reintervention Death
0 0 0 1
Beginner tip: Rows are ‘from’, columns are ‘to’. Each row must sum to 1. ‘Death’ is an absorbing row [0,0,0,1].
2.2.3 2.2 Matrix Construction for CABG
Code
# CABG Transition Matrixtrans_cabg <-matrix(c( cabg_trans$stable_to_stable, cabg_trans$stable_to_mi, cabg_trans$stable_to_reintervention, cabg_trans$stable_to_death, # From Stable... cabg_trans$mi_to_stable, 0, cabg_trans$mi_to_reintervention, cabg_trans$mi_to_death, # From MI... cabg_trans$reintervention_to_stable, 0, 0, cabg_trans$reintervention_to_death, # From Reintervention...0, 0, 0, 1# From Death (absorbing state) ),nrow =4,byrow =TRUE,dimnames =list(c("Stable", "MI", "Reintervention", "Death"), c("Stable", "MI", "Reintervention", "Death")))# Print the matrices to see their structureprint("CABG Transition Matrix:")
[1] "CABG Transition Matrix:"
Code
print(trans_cabg)
Stable MI Reintervention Death
Stable 0.8 0.08 0.07 0.05
MI 0.6 0.00 0.10 0.30
Reintervention 0.8 0.00 0.00 0.20
Death 0.0 0.00 0.00 1.00
2.2.4 2.2.1 Validation of Transition Matrices
Code
rowSums(trans_cabg)
Stable MI Reintervention Death
1 1 1 1
Code
trans_cabg["Death", ]
Stable MI Reintervention Death
0 0 0 1
Beginner tip: Rows are ‘from’, columns are ‘to’. Each row must sum to 1. ‘Death’ is an absorbing row [0,0,0,1].
2.3 Step 3: Run the Markov Model and Obtain Results
This is the core of our analysis. We will simulate a cohort of 1,000 patients over a 10-year period. We use a for loop to run the simulation for each year.
The key calculation in this step is matrix multiplication, denoted by %*% in R. In each cycle, we multiply the number of patients in each health state by the transition matrix to determine how many patients move to each new state.
2.3.1 3.1 Simulate Patient Distribution
2.3.1.1 3.1.1 Set up the simulation
We are starting with a cohort of 1,000 patients in the ‘Stable’ state.
Code
initial_pop <-c(Stable =1000, MI =0, Reintervention =0, Death =0)num_years <-10# 10 years + initial state
2.3.1.2 3.1.2 Storage for results
Code
# Storage matrices for distribution over timepci_dist <-matrix(NA_real_, nrow = num_years, ncol =length(state_names),dimnames =list(Year =1:num_years, State = state_names))cabg_dist <- pci_dist# Initial distributionpci_dist[1, ] <- initial_popcabg_dist[1, ] <- initial_pop
PCI Storage Matrix
Code
pci_dist
State
Year Stable MI Reintervention Death
1 1000 0 0 0
2 NA NA NA NA
3 NA NA NA NA
4 NA NA NA NA
5 NA NA NA NA
6 NA NA NA NA
7 NA NA NA NA
8 NA NA NA NA
9 NA NA NA NA
10 NA NA NA NA
CABG Storage Matrix
Code
cabg_dist
State
Year Stable MI Reintervention Death
1 1000 0 0 0
2 NA NA NA NA
3 NA NA NA NA
4 NA NA NA NA
5 NA NA NA NA
6 NA NA NA NA
7 NA NA NA NA
8 NA NA NA NA
9 NA NA NA NA
10 NA NA NA NA
2.3.1.3 3.1.3 Loop through years - Run the simulation and calculate state distributions
In R, a for loop is used to repeat a block of code multiple times.
Here, the loop runs once for each year of the simulation, starting from the second year (t = 2) until the final year (t = num_years).
We start at 2 because the first year (t = 1) already contains the initial patient distribution.
2.3.1.5 Step-by-Step Explanation
for (t in 2:num_years)
This sets up the loop.
On the first run, t = 2. On the next run, t = 3, and so on, until t = num_years.
pci_dist[t-1, ]
This extracts the patient distribution from the previous year (t-1) for PCI.
It is a vector showing how many patients were in each health state (Stable, MI, Reintervention, Death).
%*% trans_pci
%*% means matrix multiplication in R.
Multiplying the previous year’s patient distribution by the PCI transition matrix redistributes patients into new states for the current year.
pci_dist[t, ] <- ...
The result of the multiplication (the new patient distribution) is stored in row t of the PCI distribution matrix.
cabg_dist[t, ] <- cabg_dist[t-1, ] %*% trans_cabg
The same operation is done for CABG, using its own transition matrix.
2.3.1.6 Intuition
Each iteration of the loop moves the entire patient cohort forward by one cycle (one year).
Patients from each state “flow” into other states according to the transition probabilities.
Over the course of the loop, row by row, we build up the full history of how patients evolve under PCI and CABG.
Formula to remember:
Current year’s distribution = Previous year’s distribution × Transition matrix
2.3.1.7 Show year wise patient distribution (Trace)
This final step calculates the key metrics for our HTA: the ICER and the NMB. These metrics help us interpret the results of our model and make a decision about the cost-effectiveness of a new technology.
2.4.1 ICER (Incremental Cost-Effectiveness Ratio)
The ICER tells us the additional cost required to gain one additional unit of health (one QALY) from a new technology compared to the existing standard of care.
# Create a data frame for ICER and NMB resultsicer_nmb_df <-data.frame(Metric =c("ICER (INR per QALY)", "NMB PCI (INR)", "NMB CABG (INR)"),Value =c(icer, nmb_pci, nmb_cabg))icer_nmb_df %>%gt() %>%tab_header(title ="ICER and NMB Results") %>%fmt_currency(columns = Value,currency ="INR",decimals =2 )
ICER and NMB Results
Metric
Value
ICER (INR per QALY)
₹266,916.53
NMB PCI (INR)
₹199,890,022.52
NMB CABG (INR)
₹151,969,001.27
Code
icer_nmb_df
Metric Value
1 ICER (INR per QALY) 266916.5
2 NMB PCI (INR) 199890022.5
3 NMB CABG (INR) 151969001.3
2.5 Conclusion
In this detailed guide, we’ve built a Markov model from the ground up to compare PCI and CABG for coronary artery disease. We’ve defined health states, transition probabilities, costs, and utilities, and then simulated patient outcomes over a 10-year period. Finally, we calculated key economic evaluation metrics like the ICER and NMB to inform decision-making. This model provides a framework for understanding the long-term implications of different treatment strategies and can be adapted for other diseases and interventions. Remember, the accuracy of a Markov model depends heavily on the quality of the input data, so always ensure your parameters are based on robust clinical evidence. Happy modeling!